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ABSTRACT 

We have developed a CHi-square cOde for parameteRized modeling and char- 
acterization of photometry and Spectrophotometry (CHORIZOS). CHORIZOS 
can use up to two intrinsic free parameters (e.g. temperature and gravity for stars; 
type and redshift for galaxies; or age and metallicity for stellar clusters) and two 
extrinsic ones (amount and type of extinction). The code uses \ 2 minimization 
to find all models compatible with the observed data in the model A-dimensional 
(N = 1, 2, 3, 4) parameter space. CHORIZOS can use either correlated or uncor- 
rected colors as input and is especially designed to identify possible parameter 
degeneracies and multiple solutions. The code is written in IDL and is available 
to the astronomical community. Here we present the techniques used, test the 
code, apply it to a few well-known astronomical problems, and suggest possi- 
ble applications. As a first scientific result from CHORIZOS, we confirm from 
photometry the need for a revised temperature-spectral type scale for OB stars 
previously derived from spectroscopy. 

Subject headings: galaxies: star clusters — methods: data analysis — methods: 
numerical — methods: statistical — stars: fundamental parameters — tech- 
niques: photometric 
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1. Introduction 

A general problem in astronomy is that of finding the correspondence between the ob- 
served spectrophotometric and/or photometric properties and a series of models parameter- 
ized in terms of physical quantities. Perhaps the best well-known example is the utilization 
of Johnson U — B and B — V colors to measure the temperature and extinction of main- 
sequence stars (see, e.g. Binney & Merrifield 1998) but the problem appears in a variety of 
contexts, such as the calculation of photometric redshifts with optical/IR photometry (Koo 
1999; Benitez 2000), of extinction laws using a combination of UV-to-IR spectroscopy and 
photometry (Cardelli et al. 1989; Fitzpatrick 1999), or of stellar cluster ages and metallicities 
using broad-band colors (Girardi 2000; Whitmore et al. 1999; de Grijs et al. 2003). Each one 
of those cases has its own specific peculiarities, but they can all be considered as examples 
of the following general problem. A family of spectral energy distribution (SED) models 
f(\',Pi,P2 ■ ■ ■ ,Pn) is generated as a function of N parameters. Some of those parameters pi 
may depend on the nature and distance to the objects (e.g. temperature, age, metallicity, 
or redshift) while others depend on the properties of the intervening ISM (amount and type 
of extinction). The first type of parameters will be called in this paper intrinsic (to the 
object) and the second type extrinsic (to the object). Our data consists of one or several 
objects with measured magnitudes mi,m,2 ■ ■ - ttim+i (M > N), each with its independent 
uncertainty s±, s 2 . . . sm+i from which we can derive M independent colors 1 Ci, c 2 . . . cm- The 
general problem can be then expressed as finding what model SEDs are compatible with the 
observed colors 2 . The solution could be: 

• Unique, if the observed values for the colors and their uncertainties determine a single 
set of connected SED models (i.e. a connected iV— volume in N— dimensional parame- 
ter space) that can be described by a single set of parameters with their corresponding 
uncertainties. 

• Multiple, if the provided colors do not allow to differentiate between two or more sets 
of connected SED models. 

• Non-existent, if the properties of the object fall outside the parameter range of the 
SED models, the chosen SED models are an incorrect description of the object, or the 
photometry has systematic errors. 



1 Alternatively, we can have measured colors with their corresponding uncertainties. 

2 The photometric redshift case requires a slight reformulation of the problem because one of the param- 
eters, redshift, changes not only colors but also magnitudes. 
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Several related codes are discussed in the literature (see e.g. Romaniello et al. 2002; 
Bemtez 2000; de Grijs et al. 2003) but they all have one shortcoming: they are designed 
to deal only with an specific case of the general problem (stellar temperatures, photometric 
redhifts, cluster ages). Furthermore, some of the codes are discussed in an article but are 
not available to the astronomical community, thus hampering the testing of results by other 
groups. Finally, the treatment of how uncertainties in the measured quantities affect the 
derived parameters is in many cases poor or inexistent (some photometric redshift codes are 
an exception, see e.g. Bemtez 2000). 

We discuss in this paper CHORIZOS, a code that solves the general problem of identify- 
ing which models are compatible with an observed set of colors. In section 2 we analyze the 
problems that have to be dealt with and in section 3 we discuss the techniques to overcome 
them. In section 4 we apply an implementation of the algorithm to a number of cases and 
we end in section 5 with a summary. 

2. Description 

2.1. Problem definition 

We start by defining two different cases as a function of the number of colors and 
parameters, M = N and M > N. For M = N, we can establish as many equations (sets 
of model parameters that produce a given color) as unknowns (model parameters) and the 
problem can be treated as that of finding an exact solution to a (complicated) algebraic 
system of equations. Barring strict degeneracies in the system of equations, one can find 
either zero, one, or a finite number of solutions in N— dimensional parameter space depending 
on the specific topology of the N— volume defined in M— dimensional color space by all 
the possible parameter combinations. Adding the corresponding uncertainties associated 
with each color produces at least one connected iV— volume around each of the solutions 
in parameter space and can also generate new iV— volumes disconnected from any of the 
strict solutions and/or establish connections between solutions 3 . For M > N, the problem 
has more equations than unknowns and no exact solutions should be expected. However, 
approximate solutions that are compatible with the measured uncertainties can still be found. 
In this case, our goal should be to find that iV— volume of approximate solutions by using 
e.g x 2 minimization. 



3 Strictly speaking, once uncertainties are included, any solution is possible. In practice, we normally 
consider that the probability of a color having a real value many sigmas away from its measured one is zero. 
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2.2. M = N: the main-sequence B — V vs. U — B example 

Given the complexity and non-linearity of the general problem, different topologies can 
exist, leading to different solution types as a function of the measured magnitudes. In order 
to visualize them, we analyze the well-known example of using Johnson U — B and B — V 
colors to measure stellar temperature and reddening (M = N — 2). We show in the left 
panel of Fig. 1 the locations in 2-D color space of unreddened main-sequence stars using as 
model atmospheres those of Kurucz (2004) with Z = 0.0. The Kurucz atmospheres were 
extinguished using a Cardelli et al. (1989) law with #5495 =3.1 from £(4405 - 5495) = 0.0 to 
5.0 4 ; five different values of the temperature are shown. We also show in Fig. 2 a similar plot 
but with more temperature values (here some of the symbols have been omitted for clarity). 
Nine examples of measured magnitudes are shown in Fig. 2, each one of them with different 
measured magnitudes U, B, V, but with the same uncertainty on = Ob = oy = 0.026 in 
each case. The corresponding solutions (as calculated by CHORIZOS) are shown in Fig. 3 
as likelihood contour plots. The examples have been selected to reflect the different solution 
types. 

• Example 1 is the ideal observational situation: the measured colors correspond to a 
unique solution and the inclusion of uncertainties leads to a single set of connected 
solutions around it. We can see that because the corresponding shaded area in Fig. 2 
is crossed only by blue lines, which correspond to models with T = 9 250 — 50 000 K 
and this leads to a single connected region in Fig. 3. 

• Example 2 falls in a region where, due to the change in direction experienced by the 
zero-extinction color-color curve around T = 9 000 K, two different sets of connected 
solutions exist. Its shaded area in Fig. 2 is crossed by both blue and green lines and 
tracing them back to the zero-extinction case we arrive at two possible different ranges 
of values for the temperature. 

• Example 3 falls in the region in Fig. 2 to the upper right from where the zero-extinction 
color-color curve has changed direction for the second time. As a result, the region is 
crossed by lines corresponding to three temperature ranges (represented in Fig. 2 in 
blue, green and red, respectively) and three different sets of connected solutions are 
possible. This translates into three peaks in Fig. 3. 



4 _E(4405 — 5495) and .R5495 are the monochromatic equivalents to E(B — V) and Ry, respectively. 4405 
and 5495 are the assumed central wavelengths (in A) of the B and V filters, respectively. Monochromatic 
quantities are used because E(B — V) and Ry depend not only on the amount and type of dust but also on 
the stellar atmospheres. 
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• Examples 4, 5, and 6 correspond to the situation where the measured U — B and B — V 
colors are incompatible with any of the models but the inclusion of the corresponding 
uncertainties generates a single connected region in Fig. 3. In the case of examples 
4 and 5, the nearest edge of the N— volume of allowed colors in Fig. 2 corresponds 
to an extreme in N— dimensional parameter space (maximum temperature - 50 000 
K - for example 4, minimum £(4405 — 5495) - 0.0 - for example 5). Therefore, the 
corresponding likelihood contour plots in Fig. 3 show abrupt edges. On the other 
hand, the nearest edge of the N— volume of allowed colors for example 6 does not 
correspond to such an extreme; instead, it is caused by the change in direction of the 
zero-extinction color-color curve around T = 9 000 K. Therefore, its likelihood contour 
plot does not show an abrupt edge. 

• For example 7 the measured U — B and B — V values fall outside the M— volume of 
allowed colors but in this case there are two nearby boundaries, one of the same type 
as that of example 5 and another one of the same type as that of example 6. This 
translates into two peaks in Fig. 3, one with an abrupt edge and another one without 
it. 

• Examples 8 and 9 correspond to the cases where at least one solution exists for the 
measured colors (one for example 8, two for example 9). Here, the inclusion of uncer- 
tainties not only generates a connected region around each one of them but also a new 
one, leading to the existence of two peaks in Fig. 3 for example 8 and of three peaks 
for example 9. 

• Finally, though we have not plotted them, we can describe two possible additional 
situations. One would be the case where the shaded region in Fig. 2 was outside the 
M— volume of allowed colors and far from an edge. That would indicate that the data 
and the models are incompatible. Another situation would be reached by increasing 
the uncertainties in e.g. example 2. Then, we could have two possible solutions for 
the measured values of U — B and B — V but only a single set of connected solutions, 
since the shaded region in Fig. 2 would extend to the line defined by the turnaround 
point in temperature around 9 000 K. 

We can group our examples as a function of their solution into the three categories 
described in the introduction: unique, multiple, and non-existent. It is important to dif- 
ferentiate between cases with unique and multiple solutions for the following reason. For 
unique solutions (e.g. example 1), the derived parameters can be correctly characterized by 
a single set of mean values, T and £(4405 — 5495), and a covariance matrix defined by a T , 
0\e(44O5-5495) > an d o"t,_b(4405-5495) • 111 principle, one could do the same for multiple solutions 
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(e.g. example 2), but that would not be a correct characterization, since in those cases there 
are two or more peaks in the likelihood contour plot. In those cases, given the strong devi- 
ations from a Gaussian distribution, it would be misleading to give mean parameter values, 
since those may fall in regions of very low probability while the main peaks may be located 
at distances around or above 1 a from them (see. e.g. example 2, where the mean values 
are T = 9 910 K, £(4405 - 5495) = 0.383, and the standard deviations a T = 1 350 K and 

®E (4405 -5495) = 0.139. 

It is also important to note that, even in those cases where a unique solution is found, 
there is typically a strong correlation between the parameters. This is seen in Fig. 3 in 
that, if we were to approximate each peak by an ellipsoid, the two principal axes would 
be inclined with respect to the x and y axes. This effect is commonly referred to in the 
literature as a degeneracy between the two parameters. E.g., for example 6 we may say that 
T and £(4405 — 5495) are degenerate because our colors are compatible with a single set 
of connected solutions where likely deviations from the mean require either (a) an increase 
in both temperature and reddening or (b) a decrease in both, but not e.g. an increase in 
temperature and a decrease in reddening. 

2.3. M > N and solution existence 

The introduction of additional colors and parameters beyond two is straightforward if we 
keep M = N (though not as easy to plot). Instead of dealing with the intersection between 
ellipses and regions of a plane, we have to find the intersection between M— ellipsoids and 
regions of an N— (or M— ) dimensional space, which does not introduce any new behavior 
from the topological point of view in our description. The situation is different if M > N. 
Take as an example M = 3 and N = 2. There, we have that the two available parameters 
generate a surface (N— volume) of possible solutions inside the volume (M— volume) of all 
possible color combinations. Given the difference in dimensions, a given set of three mea- 
sured colors will always fall outside the solution surface. However, adding uncertainties to 
the measured colors will generate a finite ellipsoid (M— ellipsoid) in color space that, barring 
problems with the data or the models, should intercept the solution surface. This is repre- 
sented on the right panel of Fig. 4 for the simple case where the solution surface is a plane 
perpendicular to the major axis of the uncertainty ellipsoid. For more complex solution 
surfaces, such as the one in the left panel of Fig. 4, the intersections can have more complex 
shapes but the topological classification into unique (one connected intersection surface), 
multiple (two or more intersection surfaces), and non-existent (no intersection) solution sets 
remains unchanged. 
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Therefore, for M > N we cannot have strict solutions but only approximate ones. We 
can characterize those by defining x 2 f° r the case of uncorrelated uncertainties as: 



where c mjmo d are the model colors and a m the measured color uncertainties. For Gaussian 
uncertainties, the likelihood is then given by C = exp(— % 2 /2) and maximizing it is equivalent 
to finding the model(s) that minimize(s) x 2 ■ It can be shown (see e.g. Gould 2003) that 
the expected value of x m m is M — iV, with a standard deviation of a/2(M — N). Note that 
those are the results expected even for the case M = N, where we should find a model 
with the same exact colors as those measured and, therefore, have a Xmm °f exactly zero. 
Our definition of solution existence for M > N would be to have a Xmm reasonably close to 
M - N in units of y/2(M - N). 

An exception to the last statement should be noted. If the measured value lies close to 
an (N — 1)— dimensional edge of the solution N— volume, then the value of Xmm could be 
larger than M — N. An easy way to see this is with the most extreme case of M = N. In 
examples 4 to 7 of the previous section we saw that it is reasonable to have measured colors 
incompatible with the model ones if they fall just outside the edge of the region spanned by 
the models. In those situations where the data is close to such an edge, it is then normal 
to find somewhat higher values of x m m : f° r the specific case M = N, one should then only 
reject the existence of a solution if x m m ^ 1- 

It is still possible to find data for which there are no solutions, either for the M = N 
case or for the M > N ones. The following is a checklist of possible causes, some general 
and some specific to astronomical photometry: 

• Model validity: the models may be an incorrect description of the real SED. 

• Model applicability: the type of object or the parameter range selected may be incor- 



• Model approximations: make sure that you are deriving your colors the right way. 
Common errors are using the wrong extinction law and neglecting curvature effects in 
extinction trajectories for large values of £(4405 — 5495) (see Fig. 1). 

• Photometric quality: were the data properly acquired and calibrated? 

• Zero-point calibration: some filters require small magnitude corrections before their 
magnitudes can be used for spectrophotometry (see e.g. Cohen et al. 2003). 




(1) 



rect. 
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• Filter-throughput errors: make sure that the throughputs of your filter are correctly 
characterized. Be on the alert for filters with red/blue leaks and long tails. 

• Color transformations: beware of transforming from one filter system to another. 
Whenever possible, use the throughput definitions for the filters with which the data 
was acquired. 

3. Techniques 

Apparently, two different techniques should be needed to solve the two different cases 
defined in the previous section. In practice, both can be solved using a x 2 minimization 
algorithm. That is so because for M — N, the algebraic solution to the system of equations 
is also the solution that makes x 2 — 0, which is obviously the minimum possible value for x 2 - 
Furthermore, since £ = exp(— x 2 /2) independently of the number of degrees of freedom, we 
can use that definition to estimate the uncertainties in our parameters not only for M > N 
but also for M = N. 

CHORIZOS works by computing x 2 over the N— dimensional parameter grid and cal- 
culating the corresponding likelihood at each grid point. The current version of the code is 
written in IDL and handles up to N = 4: two parameters from the intrinsic properties of the 
SED family plus two extinction-related parameters, E(A405 — 5495) and extinction law type. 
The latter includes the i?5495-dependent family of extinction laws of Cardelli et al. (1989), 
the average LMC and LMC2 laws of Misselt et al. (1999), and the SMC law of Gordon & 
Clayton (1998). CHORIZOS starts by reading the unreddened SED models, extinguishing 
them, and obtaining the synthetic photometry at each point in a (coarse) 4-D grid. This is 
done in order to correctly deal with non-linear extinction effects. This preliminary step needs 
to be executed only once and the result can be stored in the form of binary FITS tables for 
later use. Currently, CHORIZOS includes pre-calculated tables for Kurucz (Kurucz 2004), 
Lejeune (Lejeune et al. 1997), and TLUSTY (Lanz & Hubeny 2003) stellar models and Star- 
burst 99 (Leitherer et al. 1999) cluster models using a total of 78 filter passbands (including 
the Johnson, Stromgren, 2MASS, and several HST instrument systems). The two intrinsic 
parameters are temperature and gravity for the stellar models and age and metallicity for 
the cluster models. More models and passbands will be included in the future and the user 
will also be able to add his/her own. 

The program first reads the photometry from a user-provided table, as well as a series 
of input parameters, such as which model family should be used and how fine a parameter 
grid should be employed. At this point, the user can restrict any of the four parameters 
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for any of the objects (stars, clusters or galaxies) in the input table, using either a specific 
value or a range between a minimum and a maximum. CHORIZOS then reads the model 
photometry from the previously calculated tables and interpolates to the user-selected (fine) 
grid. For each object, the likelihood is calculated at each grid point, the average values of 
the four parameters and their (output) 4x4 covariance matrix are computed, and the result 
is written in an individual file for each star. This file is read by a second program which 
produces the graphical output and allows for further manipulation of the results, such as 
calculation of stellar absolute magnitudes or cluster stellar masses. 

There are several issues regarding algorithm details and results interpretation for a 
program of this type that need to be discussed. First, if one measures magnitudes directly 
and then calculates two different colors that include one filter in common (e.g. B in the 
example discussed in the previous section), then the probability distributions for the two 
colors will be correlated (see, e.g. the data plotted in Fig. 2). It can be easily shown that 
for the case of the two colors X — Y and Y — Z formed from the three filters X, Y, and Z, 
the input (or color) covariance matrix is: 

C X -y,y-z = cov(X - Y, Y - Z) = ( ^ + / Y ~f ) , (2) 

where ax, cry, and az are the uncertainties which correspond to each filter. For correlated 
uncertainties, Eqn. 1 is no longer valid and one has to use (see e.g. Gould 2003): 

M M 

X 2 = ^2^2 i c l ~ Q,mod) &lm {Cm - C m , mo d) , (3) 
1=1 m=l 

where B = the inverse of the M x M input (or color) covariance matrix. Note, however, 
that in some cases two colors such as U — B and B — V can be considered to be in a 
first approximation as uncorrelated. That would be the case if those colors are built from 
the combination of a large number of independent measurements of U — B and B — V, a 
situation that one may encounter when collecting data from the literature. CHORIZOS 
allows either option to be used: if individual magnitudes are inputted, then the full input 
covariance matrix is utilized; if colors are chosen, then only the corresponding diagonal terms 
are used. Given the relative complexity of the algorithm implementation, We conducted an 
independent test using a Montecarlo simulation in which the measured colors were varied 
according to the input covariance matrix and the solution that minimized x 2 was selected 
in each case. The Montecarlo simulation was run with 1000 samples and the results were 
combined to produce a likelihood map which was then compared to the CHORIZOS output 
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using a variety of input data and models 5 . Results agreed in all cases. 

Another issue is that of grid size. When calculating the solution, CHORIZOS also 
checks whether the width of the resulting distribution in parameter space is comparable to 
that of the grid size near the optimum solution and warns the user if it finds that condition 
is not met so that a new run with a finer grid can be executed. A related issue is that of 
the value of xLm f° r the case M — N. As we have mentioned, it should be zero, unless 
the measured colors are incompatible with the models (e.g. the solution falls outside the 
N— volume of allowed colors). Therefore, the value of xLm can be used to check such a 
circumstance. Since we are using a grid to evaluate x 2 , the minimum value will not be zero 
but if the grid is fine enough it should be <C 1. But what if the grid is not fine enough? To 
solve this problem, after CHORIZOS finds the solution that minimizes x 2 , h interpolates 
the colors from the adjacent points into an ultrafine grid in order to obtain a more accurate 
value of the minimum. This value can then be used to decide whether the measured colors 
fall outside the N— volume of allowed colors or not. 

A problem also arises when the measured colors fall exactly at the edge of the N- volume 
of allowed colors. One example would be if in the temperature-reddening case using U — B 
and B — V the measured colors were to correspond to those of a star with zero reddening (the 
resulting likelihood diagram could be similar to that of example 5 in Fig. 3). In that case, 
the resulting mean reddening value from the likelihood data would be greater than zero and 
if we were to measure a number of stars with real i?(4405 — 5495) equal to zero, we would 
obtain positive values in all cases, leading us to an incorrect estimation of the value of the 
reddening. One solution which is implemented in CHORIZOS as an option is to extrapolate 
the grid to values beyond those provided by the original models. Using that option, the 
likelihood plot for example 5 in Fig. 3 would show the full ellipsoid and the correct mean 
value could be obtained (which in that specific case would be negative, given the location 
of point 5 in Fig. 2). This option should be used with caution, though, since extending the 
grid too much can easily lead to the inadvertent introduction of false additional solutions. 
Also, in some cases the use of extrapolated values is not recommended due to the existence 
of color near-degeneracies at the grid edges (e.g. the optical colors at the high-temperature 
end for O stars). 

Finally, the validity of the derived mean values and covariance matrix when the cal- 
culated parameter distribution is far from being a Gaussian (e.g. when multiple solutions 
exist) has to be analyzed. We will do this in the following section, where we discuss some 



5 Note that a Montecarlo method is a valid alternative way of attacking this problem but it is more 
expensive from the computational point of view. 
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sample applications of CHORIZOS. 

4. Sample applications 

4.1. Multiple solutions for stars using Johnson-Cousins photometry 

As we have seen, U—B and B — V colors alone are not enough to provide a single solution 
under all circumstances for main-sequence stellar atmospheres using a fixed extinction law. 
But what if we use additional filters in the Johnson- Cousins set? If we add the / filter and 
use V — I as a third color, we have the situation represented on the left panel in Fig. 4. We 
can see that in this case a third color eliminates most or all the multiple solutions. However, 
since we now have M = 3 > N = 2, the problem does not have an exact solution for an 
arbitrary combination of colors and we are forced into using a x 2 minimization or a similar 
technique. 

We analyze here a similar case using five Johnson- Cousins filters and a total of four 
colors: U — B, B — V, V — R, and V — I to measure the temperature and reddening of 
a main-sequence star. We use Z = 0.0, main-sequence Kurucz atmospheres, restrict the 
extinction law to be that of Cardelli et al. (1989) with -R5495 = 3.1, and fix the uncertainties 
in all five passbands to be 0.01 magnitudes. The grid was extrapolated in £(4405 — 5495) but 
not in T. In order to compare the results obtained with the full UBVRI photometry with 
those of only UBV photometry, We ran CHORIZOS using its test mode. In that mode, 
CHORIZOS is fed the true model colors and is executed to see if it is able to reproduce 
them. If there is a single well-behaved solution, the measured mean parameters should be 
very similar to the input ones; if there are multiple solutions, that may not be the case. One 
way to quantify the effect is to use the mean calculated temperature and reddening, T and 
£(4405 — 5495), and their uncertainties estimated from their standard deviations, ot and 
c r E(4405-5495)> to calculate their normalized distances from the input value: 

d T = (T-T input )/<T T , (4) 



(£(4405 - 5495) - £(4405 - 5495) input )/^ (4405 _ 5495) . (5) 

For a single well-behaved solution one expects low values of a? and o\e(4405-5495) (unless 
the input uncertainties themselves are large) and, more importantly, values of \dr\ and 
I ^£(4405-5495) I < 1-0 (and typically lower than 0.5). If there are two or more solutions, the 
normalized distances can easily be larger than one because for such a distribution (formed 
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by e.g. two well-separated narrow quasi-Gaussians) the maxima can be separated from the 
mean by more than one standard deviation. We have plotted those four quantities for the 
test case with only UBV photometry in Fig. 5. 

• For T > 15 000 K, CHORIZOS detects the existent unique solutions, as evidenced by 
the relatively low values of the uncertainties and absolute values of the normalized 
distances. Two minor points can be noted: the somewhat larger values of ot around 
37000 K are caused by the near-degeneracy of the optical colors of O stars. Also, the 
green region at the right border of the normalized distance plots is caused by edge 
effects (we did not extrapolate in temperature). 

• For T < 15 000 K, CHORIZOS produces in general much larger values of the uncer- 
tainties (especially for E(4405 — 5495)) and the normalized distances are of the order 
of or larger than one (in absolute value) for most of the region. This is caused by the 
existence of multiple solutions. 

These results are simply what we expected. The interesting point is that this demon- 
strates that when CHORIZOS is run in test mode, it can be used to predict which parameter 
ranges can be measured with the existent colors and precisions and which ones cannot be 
measured. We now turn to Fig. 6, which has the same plots but for the test case with the 
full UBVRI data using the same color scale. 

• The values of the uncertainties are much lower, being < 500 K and < 0.10 for a? and 
o\e(4405-5495) > respectively, over most of the analyzed range. The only significant devia- 
tions are those for cry around T m 40 000 K, where the introduction of RI data alleviates 
the color near-degeneracy but does not eliminate it completely, and for o"e(4405_5495) at 
certain regions near the lower left of the diagram; both effects are quite minor. 

• The values of the normalized distances (in absolute value) are also much smaller, with 
only a few regions coming close to 1.0 (in some cases due to edge effects; see, however, 
the comparison below with Lejeune atmospheres)). 

• A comparison between the likelihood plots produced by CHORIZOS from UBV and 
UBVRI cases are shown in Figs. 7 and 8. The first of those two figures shows four 
examples where the introduction of RI data yields excellent solutions, in some cases 
eliminating one of the multiple solutions and in others simply reducing the extent of 
the possible range of values. The second figure shows two of the cases where \dr\ and 
I ^£(4405-5495) I are close to 1.0 even for the UBVRI case. We can see that, even in those 
"bad" cases, the introduction of RI data is useful in further constraining the values 
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of the parameters. The reasons for the high values of the normalized distances are, in 
one case, the persistence of two solutions (with a narrowing of each peak) and, in the 
other one, the reduction to a single solution with an asymmetrical distribution. 

The importance of extending the grid for some parameters can be seen in Fig. 9, where 
we plot the values of 6^(4405-5495) for two test runs identical to the ones just described with 
the only exception of not using the grid extension option. We see there how |g?e(4405-5495) | is 
significantly larger around £"(4405 - 5495) = 0.0 and £"(4405 - 5495) = 5.0. 

In order to test this application with real data, we selected a sample of main-sequence 
stars from the list of MK standard stars by Garcia (1989) and we obtained their UBVRI 
photometry from the online catalog of Mermilliod et al. (1997). The stars were selected 
to have spectral types between B and M (an example with O stars is analyzed in the next 
subsection). Results are shown in Table 1, where the reference temperature-spectral type 
conversion has been obtained from Bessell et al. (1998) for B stars and for the rest of the 
spectral types from Carroll & Ostlie (1996). When using real data for this application, we 
are not only testing the existence of multiple solutions for a given color combination, but 
also the accuracy of the atmosphere models. For that reason, we executed CHORIZOS using 
both Kurucz and Lejeune atmospheres. 

• For A-K stars, both runs (with Kurucz and Lejeune atmospheres) produce good results, 
once the expected uncertainty in the reference temperature for a given spectral type is 
included. x mm / (M ~ N) values are low and a single solution of the correct temperature 
appears in the likelihood plots (see Fig. 10. The only significant difference between 
the two atmosphere models for these spectral types takes place for F stars, where the 
Lejeune atmospheres provide a slightly better fit than the Kurucz ones. 

• For B stars, both runs provide essentially identical results but the temperature scale 
appears to be offset by ~ 3000 K for the earliest subtypes. This could be a continuation 
of the similar effect detected by Garcia & Bianchi (2004) and other authors for O stars 
using spectroscopic data, as described in the next subsection. 

• For M stars, the Kurucz run yields bad fits (large values of xLm/ (M — N)) and clearly 
incorrect temperatures. This indicates that the Kurucz atmospheres do not provide a 
good representation for some of the optical colors of M dwarfs, as already pointed out 
by other authors (see, e.g. Lejeune et al. 1998). The Lejeune run yields low values 
of Xmm/ (M — N) (except for HD 209290) and acceptable values of T but with large 
uncertainties. The explanation for the Lejeune results can be seen in the two lower 
plots of Fig. 10: two solutions are present, the correct one around 3 650 K and another 
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Table 1. FITMODEL temperature measurements from Johnson UBVRI photometry. 



Star 


Spectral 




Temperature 






X 2 m J(M 


-AO 






type 


Reference 


Kurucz 


Lejeune 


Kurucz 


Lejeune 


HD 


36512 


B0 V 


29800 


25989± 


678 


25990± 


679 


0.17 


0.17 


HD 


37018 


Bl V 


26600 


22704± 


721 


22707± 


721 


2.39 


2.39 


HD 


74280 


B3 V 


18750 


16593± 


368 


16595± 


369 


0.12 


0.12 


HD 


219688 


B5 V 


15350 


13821± 


219 


13830± 


215 


0.15 


0.16 


HD 


222661 


B9.5 V 


10300 


10109± 


206 


10135± 


129 


0.09 


0.02 


HD 


18331 


Al V 


9230 


9239± 


405 


9097± 


528 


0.78 


0.50 


HD 


216956 


A3 V 


8720 


8812± 


235 


8595± 


402 


1.45 


1.41 


HD 


26911 


F3 V 


6700 


7417± 


75 


6988± 


70 


0.88 


0.64 


HD 


27534 


F5 V 


6440 


7226± 


95 


6811± 


98 


1.14 


0.02 


HD 


222368 


F7 V 


6300 


6707± 


243 


6321± 


202 


0.43 


0.38 


HD 


102870 


F9 V 


6100 


5906± 


227 


5945± 


545 


2.27 


0.90 


HD 


141004 


GO V 


6030 


5895± 


76 


5799± 


368 


0.67 


0.62 


HD 


27836 


Gl V 


5950 


6566± 


622 


6524± 


418 


3.02 


0.81 


HD 


1835 


G2.5 V 


5830 


5547± 


44 


5399± 


37 


0.55 


0.32 


HD 


20630 


G5 V 


5770 


5722± 


49 


5529± 


37 


0.02 


1.03 


HD 


20794 


G8 V 


5570 


5776± 


52 


5570± 


71 


0.35 


0.30 


HD 


26965 


K0.5 V 


5160 


5363± 


36 


5263± 


26 


0.92 


0.15 


HD 


22049 


K2 V 


4900 


5153± 


34 


5122± 


21 


0.22 


0.44 


HD 


196795 


K5 V 


4350 


4335± 


454 


4021± 


42 


0.53 


2.15 


HD 


209290 


M0.5 V 


3780 


5303± 


143 


4714± 


722 


3.94 


3.27 


HD 


131976 


Ml. 5 V 


3650 


10080± 


188 


5975±3134 


3.10 


0.14 


HD 


36395 


Ml. 5 V 


3650 


10083± 


168 


5812±3070 


3.83 


0.19 


HD 


119850 


M2 V 


3580 


10880± 


566 


8924±3139 


3.09 


0.92 
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one around 10 000 K. This implies that the Lejeune atmospheres provide a better 
representation of the colors of M dwarfs but that UBVRI photometry alone is not 
enough to distinguish between M stars and reddened late-B stars. 

The conclusion is that the addition of RI data to UBV photometry is very useful 
in constraining the temperature and extinction of stars (under certain assumptions such 
as knowledge of the extinction law and the luminosity type) and that CHORIZOS can be 
successfully used to derive those parameters. A possible application of the code would be to 
automatically generate temperatures and extinctions from photometric surveys such as the 
one planned with the GAIA mission. 

4.2. Measuring optical-IR extinction laws 

As a second example, we use Johnson UBV and 2MASS JHK S to determine the ex- 
tinction and extinction law experienced by an early-type star. We select as input a 35 000 K, 
log p = 5.0, solar metallicity Kurucz atmosphere model with uncertainties of 0.01 magnitudes 
in each of the six filters. The Kurucz atmosphere is extinguished from £(4405 — 5495)= 0.0 
to £(4405 - 5495)= 5.0 using Cardelli et al. (1989) laws from i? 5 495 = 2.0 to i? 5 495 = 6.0. 

We first assume that we have an accurate spectral type for the star and, therefore, that 
we know a priori the temperature and gravity of the star. We do this by constraining the 
temperature and gravity in CHORIZOS to a single value (the true one). Grid extension is 
used for both £(4405 — 5495) and -R5495. Results for this M = 5, N = 2 case are shown in 
Fig. 11. The two lower plots demonstrate that UBV JHK S is an adequate choice of filters 
to measure the extinction and extinction law experienced by hot stars. The values of both 
normalized distances are very close to 0.0 everywhere with the only exception of the region 
around £(4405 — 5495)= for dn 5495 . The latter is an expected behavior, since for low values 
of the reddening all extinction laws produce similar results and for £(4405 — 5495)= 0.0 they 
are strictly degenerate. This is evidenced in the upper right plot of Fig. 11: the value of 
<jR 54ga is kept lower than 0.10 for £(4405 — 5495) > 0.4 but increases rapidly as we approach 
£(4405 — 5495) = 0.0, where the lack of information on the extinction law provided by the 
data manifests itself in large values of ctr 5495 . Note, however, that <t_b (4405 -5495) is n °t strongly 
affected by this, since its value is kept below 0.0075 everywhere. 

Suppose now that we know that the star is of early type but we cannot constrain its 
temperature farther than that due to the lack of an accurate spectral type. We can simulate 
such a case in CHORIZOS by leaving the temperature unconstrained and selecting main- 
sequence Kurucz models. Results for this M = 5, N = 3 case are shown in Fig. 12. We 
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see in the two lower plots there that the normalized distances are slightly worse than in the 
previous case, but still within acceptable ranges. The degeneracy in -R5495 is still present for 
£(4405 — 5495)= 0.0, as expected, but the information in the photometry is accurate enough 
to yield good estimates of both £(4405 — 5495) and -R5495. The loss of information caused 
by the unconstrained temperature only translates into larger uncertainties in the measured 
quantities, but even those are smaller than 0.03 anywhere for <7£(4405-5495) and lower than 
0.25 for a R5495 for £(4405 - 5495) > 0.4. 

We also tested the measurement of optical-IR extinction laws with CHORIZOS using 
real data from the literature. We selected from the sample used by Cardelli et al. (1989) 
to derive their extinction law the four stars present in the Galactic O star catalog of Mafz- 
Apellaniz et al. (2004) which have (a) UBVJHK S data in the catalog and (b) values of 
E(B - V) measured by Cardelli et al. (1989) greater than 0.45. CHORIZOS was run using 
the UBVJHKs photometry (M = 5) and TLUSTY atmospheres twice, first constraining 
the temperatures and gravities to fixed values (N = 2) and then leaving them unconstrained 
(N = 4). The values for the temperatures and gravities were derived from the spectral types 
by using a scale intermediate between the ones proposed by Vacca et al. (1996) and Garcia & 
Bianchi (2004) and then selecting the closest TLUSTY model. Results are shown in Table 2. 

• CHORIZOS results for -R5495 are in excellent agreement with those of Cardelli et al. 
(1989). That article does not provide error estimates but our results are always within 
one sigma and are also small enough for the output to be meaningful. 

• CHORIZOS results for £(4405 — 5495) are also in very good agreement with the ref- 
erence values, though the unconstrained results are in all cases lower than the ones 
provided by Cardelli et al. (1989) (but always within two sigmas). The likely origin of 
this minor difference is the use of values between —0.30 and —0.32 for the (B — V) 
colors of O stars by those authors. TLUSTY atmospheres predict redder colors by 
?s 0.02 magnitudes. Once that difference is included, the agreement is excellent. 

• The Xmin values indicate very good fits in all cases except for HD 73882, which is still 
acceptable. 

• Assuming that no biases are present (e.g. systematic errors in the TLUSTY atmo- 
spheric models), the unconstrained results favor the lower temperature scale of Garcia 
& Bianchi 2004 (04 dwarfs around 40 000 K and 06 dwarfs around 32 000-35 000 K, 
implying a boundary between O and B dwarfs below 30 000 K) over that of Vacca et al. 
1996 (with 04 dwarfs close to 50 000 K and a boundary between O and B dwarfs at 
34 000 K). 



Table 2. Comparison between FITMODEL and Cardelli et al. (1989) results for the reddening and extinction law of 

4 O stars. 



Star 


Spectral 


Cardelli et al. (1989) 




Constrained T and log g (M 


= 5, N = 2) 




Unconstrained T and logg (M = 


5, N = 4) 






type 


E(B - V) 


R v 


T (K) 


logs 


£(4405 - 5495) 


-R5495 


X 2 ■ 

M-N 


T (K) 


logs 


E(4405 - 5495) 


-R5495 


X 2 ■ 

A min 

M-N 


HD 46202 


09 V 


0.47 


3.12 


32 500 


4.00 


0.478±0.013 


3.05±0.12 


1.70 


28 900±1800 


4.13±0.46 


0.448±0.020 


3.14±0.14 


0.57 


HD 73882 


08.5 V ((n)) 


0.72 


3.39 


32 500 


4.00 


0.704±0.017 


3.42±0.10 


3.78 


28 600±1900 


4.21±0.43 


0.679±0.021 


3.44±0.10 


5.64 


HD 229196 


06 III (n)(f) 


1.22 


3.12 


37 500 


3.75 


1.221±0.014 


3.16±0.05 


0.47 


37 400±3 400 


3.88±0.56 


1.214±0.028 


3.16±0.06 


1.32 


CPD -59 2600 


06 V ((f)) 


0.53 


4.17 


37 500 


4.00 


0.507±0.023 


4.13±0.21 


1.14 


36 500±3 500 


3.90±0.56 


0.492±0.033 


4.14±0.23 


2.17 



I 

I— ' 

I 
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In summary, CHORIZOS can be used to measure reddenings and extinction laws with 
good precision, even when accurate spectral types are not available. 

4.3. What precision is required to measure gravity for O stars with optical 

photometry alone? 

For our third example we want to investigate the possibility of using Stromgren pho- 
tometry to measure the surface gravity of stars of unknown temperature and extinction 
(but with a known extinction law). This is obviously a difficult task, since O-star optical 
colors are quasi-degenerate in temperature and even more so in gravity. Our goal will be to 
determine what kind of photometric precision would be required and to assess whether such 
an accuracy is attainable. 

We use TLUSTY atmospheres with solar metallicity, T = 35 000 K, and log g between 
3.25 and 4.75 observed with Stromgren ubvy photometry. The extinction law is restricted to 
be of Cardelli et al. (1989) type with -R5495 = 3.1 but no constraints are placed on the possible 
values of T, log g, or E(4405 — 5495), yielding an M = 3, iV = 3 case. Grid extension is 
used for E(4405 — 5495) but not for T or log g due to the incompleteness of the model grid 
in these two last parameters (Lanz & Hubeny 2003). 

In a first run, values of 0.003 were used for the uncertainties in the measured magnitudes. 
Results for a\ ogg and d\ ogg are shown in the left panels of Fig. 13. We see there that the 
supplied photometry does not yield enough information to measure gravities. The high values 
of both d\ ogg and o"i ogs are characteristic of an almost constant output result of ~ 4.0 ± 0.5 
for logg independent of the input values for the gravity and reddening. 

In a second run we use values of 0.001 for the uncertainties. Results are shown in the 
right panels of Fig. 13. Now the values for d\ ogg and o"io gs show that some discrimination 
is possible among gravities, with typical uncertainties in log g around 0.10 — 0.20, which is 
enough to differentiate between main sequence stars and supergiants. 

A hypothetical observer should now ask him/herself: Am I convinced that the atmo- 
sphere models are correct to within 0.001 magnitudes? Can I calibrate the photometry to 
levels good enough to accurately measure differences of 1 milimagnitude? At the current 
level of knowledge and technology the answers to both questions are likely to be no, hence 
it can be deduced that at the present time optical photometry alone cannot be used to ac- 
curately measure O-star gravities. However, future improvements in atmosphere modeling 
and photometric accuracies may change the situation. 
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5. Summary 

CHORIZOS is a multi-purpose ^-minimization SED-fitting code that can be applied 
to either photometric or spectrophotometric data. In this article we have described the 
techniques it employs and applied it to several astronomical examples. At the time of this 
writing, a beta version of the program with a limited number of SED models and 78 filter 
passbands is available from http://www.stsci.edu/~jmaiz. In the future, a full version 
will be made available which will include user-defined filter sets or wavelength ranges as well 
as the possibility of adding other SED models. 

Besides the obvious application of selecting the model (s) which are compatible with a 
given set of observed data, CHORIZOS can be used for a number of other astronomical 
applications. It can be utilized to select the optimum choice of filters and minimum S/N 
requires when planning an observation, to test atmospheric models and extinction laws, or 
to calibrate the zero points of a filter system, to name a few. 

I would like to thank Leonardo Ubeda for his help with the testing of the code, an anony- 
mous referee for his/her helpful comments, and Rodolfo Xeneize Barba for his christening 
suggestions. 

REFERENCES 

Benitez, N. 2000, ApJ, 536, 571 

Bessell, M. S., Castelli, F., & Plez, B. 1998, A&A, 333, 231 

Binney, J., & Merrifield, M. 1998, Galactic Astronomy (Princeton: Princeton University 
Press) 

Cardelli, J. A., Clayton, G. C, & Mathis, J. S. 1989, ApJ, 345, 245 

Carroll, B. W., & Ostlie, D. A. 1996, An introduction to modern astrophysics (Addison- 
Wesley) 

Cohen, M., Wheaton, W. A., & Megeath, S. T. 2003, AJ, 126, 1090 

de Grijs, R., Fritze-v. Alvensleben, U., Anders, P., Gallagher, J. S., Bastian, N., Taylor, 
V. A., & Windhorst, R. A. 2003, MNRAS, 342, 259 

Fitzpatrick, E. L. 1999, PASP, 111, 63 



-20- 



U-B 



U-B 





• 


unred. models AT=5000 K 




A 


AE(4405-5495)-0.1 , T=50000 K 




■ 


AE(4405-5495)-0.1 , T=20000 K 




+ 


AE(4405-5495)-0.1 , T=12000 K 




★ 


AE(4405-5495)-0.1 , T=10000 K 




* 


AE(4405-5495)-0.1 , T= 5000 K 





- 3 



in 
i 

, < 



• 


unred. models AT=5000 K 








▲ 


AE(4405-5495)=0.1, R 


=2.0 




1 


■ 


AE(4405-5495)=0.1, R 


=3.0 






+ 


AE(4405-5495)=0.1 , R 


=4.0 






★ 


AE (4405-5495)=0.1, R M9 


=5.0 






# 


AE (4405-5495)=0.1, R M9 


=6.0 








-1.0 



-0.5 



0.0 



0.5 



1.0 



U-B 



U-B 



Fig. 1. — (left) U — B vs. B — V color-color plot for Z = 0.0, main- sequence Kurucz 
atmospheric models. The line with circles indicates the location of the unredenned values as 
a function of temperature, starting at T = 50 000 K, with the circles marking those points 
where the temperature is a multiple of 5 000 K. The rest of the lines indicate the colors as a 
function of reddening using the Cardelli et al. (1989) law with -R5495 = 3.1 for five different 
temperatures. Symbols are plotted at intervals of Ai?(4405 — 5495) = 0.10. (right) Same 
as left but now plotting colors as a function of reddening for five different values of -R5495 
using the same family of extinction laws for a temperature of 50 000 K. The scale has been 
enlarged to include all unredenned models with T = 3 500 — 50 000 K and all 50 000 K models 
up to £(4405 - 5495) = 5.0. 
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Fig. 2. — 5 — V vs. U — B color-color plot for the same conditions as in Fig. 1. Here we 
include a higher number of extinguished models using a color code to differentiate among 
temperature ranges which are relevant to determine the number of possible temperature + 
reddening solutions for a given (U — B) + (B — V) color combination. Nine examples are 
marked, each one of them with oy = ob = oy = 0.026. 
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Fig. 3. — Likelihood contour plots produced by CHORIZOS for the nine examples described 
in the text and shown in Fig. 2. The x axis corresponds to the temperature in K and the y 
axis to £(4405 - 5495). 
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Fig. 4. — (left) V — I vs. B — V vs. U — B 3-color plot for Z = 0.0, main-sequence Kurucz 
atmospheric models. The range plotted covers T = 3 500 — 50 000 and £"(4405 — 5495) = 
0.0 — 2.0 and the extinction law used is that of Cardelli et al. (1989) with -R5495 = 3.1. 
The color surface marks the location in 3-color space while the black one is the projection 
onto the (U — B)-(B — V) plane (see Figs. 1 and 2). (right) Basic topology for an M = 3, 
N = 2 case such as the one on the left panel. Given that M > N, the measured colors (blue 
circle) always lies outside the solution surface (red grid). However, including the uncertainty 
ellipsoid (blue grid), yields an intersection surface (solid yellow) of likely solutions. 
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Fig. 5. — (top) Plots of <jt (left) and o\e(4405-5495) (right) for sample application 1 (main- 
sequence Kurucz models with -R5495 = 3.1) as a function of logT (in K) and _E(4405 — 5495) 
calculated for UBV photometry with uncertainties of 0.01 magnitudes for each filter, (bot- 
tom) Plots of dx (left) and g?£(4405-5495) (right) corresponding to the same case. 
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Same as Fig. 5 for UBVRI photometry. 
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Fig. 7. — Likelihood contour diagrams for four selected cases in sample application 1 using 
UBV data (left) and UBVRI data (right). Temperature (in K) is plotted in the horizontal 
axis and i?(4405 — 5495) in the vertical axis. A white star is used to indicate the true 
temperature and reddening. The irregularities in the contour diagram for the third star are 
caused by grid-size effects. 




Fig. 9. — Plots of <iE(4405_5495) for sample application 1 without grid extrapolation in 
£(4405 - 5495) for the UBV (left) and UBVRI (right) cases. 
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Fig. 10. — Likelihood contour diagrams for four of the stars in Table 1 using Lejeune at- 
mospheres. From left to right and top to bottom, the cases shown are HD 18331 (Al V), 
HD 20794 (G8 V), HD 131976 (Ml. 5 V), and HD 36395 (M1.5 V). Temperature (in K) is 
plotted in the horizontal axis and i?(4405 — 5495) in the vertical axis. 
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Fig. 11. — (top) Plots of cte(4405_5495) (left) and (Jr 5495 (right) for sample application 2 (main- 
sequence 35 000 K Kurucz model) as a function of i?(4405 — 5495) and -R5495 calculated for 
UBVJHK photometry with uncertainties of 0.01 magnitudes for each filter and restricting 
the model temperature to its true value, (bottom) Plots of 6^(4405-5495) (left) and 6?# 5495 
(right) corresponding to the same case. 
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g. 12. — Same as Fig. 11 but without placing constraints on the model temperature. 
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Fig. 13. — (top) Plots of o"i og9 for sample application 3 (TLUSTY solar metallicity 35 000 
K models) as a function of log g and i^(4405 — 5495) calculated for Stromgren photometry 
with uncertainties of 0.003 (left) and 0.001 (right) magnitudes for each filter, (bottom) Plots 
of d\ ogg corresponding to the same cases. 
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